pcp_claims <- fread("~/pcp_claims.csv")

#######################################################################
###############                TABLE 2                 ################
#######################################################################

#MH
mh_mh <- felm(MH_dummy_cum3_v1~FE_MH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
ph_mh <- felm(PH_dummy_cum3_v1~FE_MH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
acsc_mh <- felm(ACSC_cum3~FE_MH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

#PH
mh_ph <- felm(MH_dummy_cum3_v1~FE_PH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
ph_ph <- felm(PH_dummy_cum3_v1~FE_PH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
acsc_ph <- felm(ACSC_cum3~FE_PH_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

#ACSC
mh_acsc <- felm(MH_dummy_cum3_v1~FE_ACSC_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
ph_acsc <- felm(PH_dummy_cum3_v1~FE_ACSC_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
acsc_acsc <- felm(ACSC_cum3~FE_ACSC_std+RaceEthnicity+Medicare+Medicaid+PH_dummy_lag1_v1+MH_dummy_lag1_v1+ACSC_lag1|Age_bins+Married+EnrollmentPriority+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

stargazer(mh_mh, mh_ph, mh_acsc, ph_mh, ph_ph, ph_acsc, acsc_mh, acsc_ph, acsc_acsc,
          apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_std", "FE_PH_std", "FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
          digits=2,
          add.lines=list(c("Mean Dep. Var. (x100)", 100*round(mean(mh_mh$fitted.values),4), 100*round(mean(ph_mh$fitted.values),4), 100*round(mean(acsc_mh$fitted.values),4))))

rm(mh_mh, mh_ph, mh_acsc, ph_mh, ph_ph, ph_acsc, acsc_mh, acsc_ph, acsc_acsc)


##############################################################
##############              TABLE 3             ############## 
##############################################################

pcp_claims$herc1 <- log(1+pcp_claims$herc_yr1)
pcp_claims$herc3 <- log(1+pcp_claims$herc_cum3)

death <- felm(Death3Year~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc1 <- felm(herc1~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>12])
herc3 <- felm(herc3~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

mh <- stargazer(death, herc1, herc3, 
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("3 Year Mortality", "Log 1Y Avg Cost", "Log 3Y Avg Cost", "3 Year ACSC"))

death <- felm(Death3Year~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc1 <- felm(herc1~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>12])
herc3 <- felm(herc3~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

ph <- stargazer(death, herc1, herc3,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2)

death <- felm(Death3Year~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
herc1 <- felm(herc1~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>12])
herc3 <- felm(herc3~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

acsc <- stargazer(death, herc1, herc3,
                  apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=2)

rbind(mh[1:17], ph[15:17], acsc[15:24])
rm(mh, ph, acsc)




####################################################################
##############                TABLE 4                ############### 
####################################################################

# Note that cause of death data requires additional data clearance from VA
cod <- read_csv("~/cod.csv")
cod$scrssn <- as.numeric(cod$scrssn); pcp_claims$ScrSSN <- as.numeric(pcp_claims$ScrSSN)
pcp_claims <- merge(x=pcp_claims, y=cod, by.x="ScrSSN", by.y="scrssn", all.x=T)

# If someone has NOT died from above, then zero
# BE CAREFUL: could mean we don't have COD data yet)
pcp_claims <- data.frame(pcp_claims)
replace_index <- which(names(pcp_claims) %in% c("Overdose", "Accident", "Cancer", "Suicide", "HeartDisease", "ChronicLiverDisease", "LowerRespiratory", "Cerebrovascular"))
pcp_claims[replace_index][is.na(pcp_claims[replace_index])] <- 0
rm(replace_index)
pcp_claims <- data.table(pcp_claims)

pcp_claims$SuspectedSuicide <- pmax(pcp_claims$Suicide, pcp_claims$Overdose, pcp_claims$Accident)

pcp_claims$Overdose_3yr <- (pcp_claims$Overdose==1 & pcp_claims$Death3Year==1)
pcp_claims$Cancer_3yr <- (pcp_claims$Cancer==1 & pcp_claims$Death3Year==1)
pcp_claims$Suicide_3yr <- (pcp_claims$Suicide==1 & pcp_claims$Death3Year==1)
pcp_claims$SuspectedSuicide_3yr <- (pcp_claims$SuspectedSuicide==1 & pcp_claims$Death3Year==1)
pcp_claims$HeartDisease_3yr <- (pcp_claims$HeartDisease==1 & pcp_claims$Death3Year==1)
pcp_claims$ChronicLiverDisease_3yr <- (pcp_claims$ChronicLiverDisease==1 & pcp_claims$Death3Year==1)
pcp_claims$LowerRespiratory_3yr <- (pcp_claims$LowerRespiratory==1 & pcp_claims$Death3Year==1)
pcp_claims$Cerebrovascular_3yr <- (pcp_claims$Cerebrovascular==1 & pcp_claims$Death3Year==1)

cancer <- felm(Cancer_3yr~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
               +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
heart <- felm(HeartDisease_3yr~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
suicide <- felm(Suicide_3yr~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
possiblesuicide <- felm(SuspectedSuicide_3yr~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                        +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
lowresp <- felm(LowerRespiratory_3yr~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
cerebro <- felm(Cerebrovascular_3yr~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

mh <- stargazer(cancer, heart, suicide, possiblesuicide, lowresp, cerebro,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=3, dep.var.labels=c("Cancer", "Heart", "Suicide", "Possible Suicides", "Lower Respiratory", "Cerebrovascular"))

cancer <- felm(Cancer_3yr~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
               +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
heart <- felm(HeartDisease_3yr~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
suicide <- felm(Suicide_3yr~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
possiblesuicide <- felm(SuspectedSuicide_3yr~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                        +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
lowresp <- felm(LowerRespiratory_3yr~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
cerebro <- felm(Cerebrovascular_3yr~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

ph <- stargazer(cancer, heart, suicide, possiblesuicide, lowresp, cerebro,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=3, dep.var.labels=c("Cancer", "Heart", "Suicide", "Possible Suicides", "Lower Respiratory", "Cerebrovascular"))

cancer <- felm(Cancer_3yr~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
               +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
heart <- felm(HeartDisease_3yr~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
suicide <- felm(Suicide_3yr~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
possiblesuicide <- felm(SuspectedSuicide_3yr~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                        +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
lowresp <- felm(LowerRespiratory_3yr~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
cerebro <- felm(Cerebrovascular_3yr~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

acsc <- stargazer(cancer, heart, suicide, possiblesuicide, lowresp, cerebro,
                  apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=3, dep.var.labels=c("Cancer", "Heart", "Suicide", "Possible Suicides", "Lower Respiratory", "Cerebrovascular"))

rbind(mh[1:17], ph[15:17], acsc[15:24])
rm(mh, ph, acsc, cancer, heart, suicide, possiblesuicide, lowresp, cerebro)



##################################################################
################            TABLE 5           ####################
##################################################################

# <5mins
visit_composition <- sqlQuery(db.handle,
                              "
set nocount on;
set ansi_warnings off;

SELECT
  PatientICN, COUNT(DISTINCT CASE WHEN type='PC' THEN VisitDate ELSE NULL END) AS N_PC, COUNT(DISTINCT CASE WHEN type='MH' THEN VisitDate ELSE NULL END) AS N_MH, 
  COUNT(DISTINCT CASE WHEN type='ED' THEN VisitDate ELSE NULL END) AS N_ED, COUNT(DISTINCT CASE WHEN type='INP' THEN VisitDate ELSE NULL END) AS N_INP, COUNT(DISTINCT VisitDate) AS N_All
FROM(
SELECT DISTINCT
  cohort.PatientICN, CAST(visit.VisitDateTime AS DATE) AS VisitDate, CASE WHEN StopCode IN ('323', '322', '342', '348', '350', '702', '160', '502', '301') THEN 'PC' WHEN (StopCode>=502 AND StopCode<=599) THEN 'MH' WHEN StopCode='130' THEN 'ED' ELSE 'OTHER' END AS type
FROM CDWWork.Outpat.Visit AS visit
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.VisitDateTime)<=(365)
INNER JOIN CDWWork.Dim.StopCode AS sc
  ON visit.PrimaryStopCodeSID=sc.StopCodeSID
WHERE visit.VisitDateTime>=CAST('2005-01-01' AS datetime2(0)) AND Visit.VisitDateTime<CAST('2018-03-01' AS datetime2(0)) AND StopCodeName NOT LIKE '%ADMIN%' AND StopCodeName NOT LIKE '%TELE%'
UNION(
SELECT
  cohort.PatientICN, CAST(AdmitDateTime AS DATE) AS VisitDate, 'INP' AS type
FROM CDWWork.Inpat.Inpatient AS visit
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.AdmitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.AdmitDateTime)<=(365)
WHERE visit.AdmitDateTime>=CAST('2005-01-01' AS datetime2(0)) AND visit.AdmitDateTime<CAST('2018-03-01' AS datetime2(0))
)
) AS z
GROUP BY PatientICN")


visit_composition <- data.table(visit_composition)
pcp_claims <- merge(x=pcp_claims, y=visit_composition, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$N_PC[which(is.na(pcp_claims$N_PC))] <- 0
pcp_claims$N_PC[which(pcp_claims$N_PC>=24)] <- 24
pcp_claims$N_MH[which(is.na(pcp_claims$N_MH))] <- 0
pcp_claims$N_MH[which(pcp_claims$N_MH>=24)] <- 24
pcp_claims$N_ED[which(is.na(pcp_claims$N_ED))] <- 0
pcp_claims$N_ED[which(pcp_claims$N_ED>=24)] <- 24
pcp_claims$N_INP[which(is.na(pcp_claims$N_INP))] <- 0
pcp_claims$N_INP[which(pcp_claims$N_INP>=10)] <- 10

pcp_claims$N_All[which(is.na(pcp_claims$N_All))] <- 0
pcp_claims$N_All[which(pcp_claims$N_All>=72)] <- 72


cms <- sqlQuery(db.handle, "SELECT DISTINCT PatientICN, VisitDateTime FROM OMHO_QFR.ECON.smc_pcpvalueadd WHERE type='OUT'")

cms <- data.table(cms)
cms$VisitDateTime <- as.Date(cms$VisitDateTime)
cms <- merge(x=cms, y=pcp_claims[,c("PatientICN", "VisitDate")], by.x="PatientICN", by.y="PatientICN")
cms$VisitDate <- as.Date(cms$VisitDate)

cms <- cms[VisitDateTime>=VisitDate & VisitDate+365>=VisitDateTime]
cms <- cms[,.(CMS_Days=length(unique(VisitDateTime))), by="PatientICN"]

pcp_claims <- merge(x=pcp_claims, y=cms[,c("PatientICN", "CMS_Days")], by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$CMS_Days[which(is.na(pcp_claims$CMS_Days))] <- 0
pcp_claims$CMS_Days[which(pcp_claims$CMS_Days>=24)] <- 24

N_All <- felm(N_All~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_PC <- felm(N_PC~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
             +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_MH <- felm(N_MH~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
             +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_ED <- felm(N_ED~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
             +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_INP <- felm(N_INP~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_CMS <- felm(CMS_Days~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[Age>=65])

mh <- stargazer(N_All, N_PC, N_MH, N_ED, N_INP, N_CMS,
                keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=3, dep.var.labels=c("All", "Primary Care", "Mental Health", "Emergency", "Inpatient"))


N_All <- felm(N_All~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_PC <- felm(N_PC~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
             +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_MH <- felm(N_MH~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
             +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_ED <- felm(N_ED~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
             +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_INP <- felm(N_INP~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_CMS <- felm(CMS_Days~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[Age>=65])

ph <- stargazer(N_All, N_PC, N_MH, N_ED, N_INP, N_CMS,
                keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=3, dep.var.labels=c("All", "Primary Care", "Mental Health", "Emergency", "Inpatient"))


N_All <- felm(N_All~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_PC <- felm(N_PC~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
             +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_MH <- felm(N_MH~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
             +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_ED <- felm(N_ED~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
             +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_INP <- felm(N_INP~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_CMS <- felm(CMS_Days~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
              +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[Age>=65])

acsc <- stargazer(N_All, N_PC, N_MH, N_ED, N_INP, N_CMS,
                  keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=3, dep.var.labels=c("All", "Primary Care", "Mental Health", "Emergency", "Inpatient"))

rbind(mh[1:17], ph[15:17], acsc[15:24])
rm(N_All, N_PC, N_MH, N_ED, N_INP, N_CMS, visit_composition, cms)



############################################################################
################                  TABLE 6               ####################
############################################################################

#All referrals from any provider
# Very fast
referrals <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;

SELECT DISTINCT
  cohort.PatientICN,
  CASE WHEN con.ConsultSID IS NOT NULL THEN 1 ELSE 0 END AS Ref,
  CASE WHEN (sc.StopCode>=502 AND sc.StopCode<=599) OR (ToRequestServiceName LIKE '%MH%' OR ToRequestServiceName LIKE '%MENTAL%' OR ToRequestServiceName LIKE '%PSYCH%' OR ToRequestServiceName LIKE '%PTSD%' OR ToRequestServiceName LIKE '%SUBST%') THEN 1 ELSE 0 END AS Ref_MH,
  CASE WHEN ToRequestServiceName LIKE '%CARDIO%' OR ToRequestServiceName LIKE '%CARDIAC%' OR ToRequestServiceName LIKE '%EKG%' OR ToRequestServiceName LIKE '%ECHO%' OR StopCodeName LIKE '%CARDIO%' OR StopCodeName LIKE '%CARDIAC%' OR StopCodeName LIKE '%EKG%' OR StopCodeName LIKE '%ECHO%' THEN 1 ELSE 0 END AS Ref_PH,
  CASE WHEN ToRequestServiceName LIKE '%CANCER%' OR ToRequestServiceName LIKE '%ONCOLOGY%' OR ToRequestServiceName LIKE '%TUMOR%' OR StopCodeName LIKE '%CANCER%' OR StopCodeName LIKE '%ONCOLOGY%'THEN 1 ELSE 0 END AS Ref_Cancer
FROM CDWWork.Con.Consult AS con
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON con.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
ON cohort.PatientICN=pat.PatientICN AND con.RequestDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, con.RequestDateTime)<=(365)
INNER JOIN CDWWork.Dim.RequestService AS req
  ON con.ToRequestServiceSID=req.RequestServiceSID
LEFT OUTER JOIN CDWWork.Dim.AssociatedStopCode AS sc
  ON con.ToRequestServiceSID=sc.RequestServiceSID
INNER JOIN CDWWork.Dim.StopCode AS sc2
  ON sc.StopCodeSID=sc2.StopCodeSID
  WHERE RequestDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND RequestDateTime<CAST('2018-06-08 00:00:00' AS datetime2(0))
  AND CPRSStatus NOT IN ('*%', 'DISCONTINUED', 'CANCELLED') 
  AND RequestType='C' -- a consult and not a procedure
  AND (AdministrativeFlag IS NULL OR NOT AdministrativeFlag='Y')
  AND InpatOutpat='O'
")

referrals <- data.table(referrals)
referrals <- referrals[PatientICN %in% pcp_claims$PatientICN]
referrals <- referrals[, .(Ref=max(Ref), Ref_MH=max(Ref_MH), Ref_PH=max(Ref_PH)), by="PatientICN"]

pcp_claims <- merge(x=pcp_claims, y=referrals, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$Ref[which(is.na(pcp_claims$Ref))] <- 0
pcp_claims$Ref_MH[which(is.na(pcp_claims$Ref_MH))] <- 0
pcp_claims$Ref_PH[which(is.na(pcp_claims$Ref_PH))] <- 0

laboratory <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;

SELECT 
  cohort.PatientICN, COUNT(DISTINCT LabChemTestSID) AS N_specimen, COUNT(DISTINCT LabPanelSID) AS N_panel
FROM CDWWork.Chem.PatientLabChem AS lab
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON lab.PatientSID=pat.PatientSID
INNER JOIN CDWWork.Dim.Location AS loc
  ON lab.RequestingLocationSID=loc.LocationSID
INNER JOIN CDWWork.Dim.StopCode AS sc
  ON loc.PrimaryStopCodeSID=sc.StopCodeSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND lab.LabChemSpecimenDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, lab.LabChemSpecimenDateTime)<=(365)
WHERE lab.LabChemSpecimenDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND lab.LabChemSpecimenDateTime<CAST('2018-06-08 00:00:00' AS datetime2(0))
AND LocationType='CLINIC' AND NOT StopCode='130'
GROUP BY cohort.PatientICN
")

laboratory <- data.table(laboratory)
pcp_claims <- merge(x=pcp_claims, y=laboratory, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$N_specimen[which(is.na(pcp_claims$N_specimen))] <- 0
pcp_claims$N_specimen[which(pcp_claims$N_specimen>=100)] <- 100
pcp_claims$N_panel[which(is.na(pcp_claims$N_panel))] <- 0
pcp_claims$N_panel[which(pcp_claims$N_panel>=30)] <- 30

## Imaging in outpatient settings 
radiology <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;

SELECT 
  cohort.PatientICN, COUNT(DISTINCT RadiologyCPTProcedureCode) AS N_imaging
FROM CDWWork.Rad.RadiologyExam AS rad
INNER JOIN CDWWork.Dim.RadiologyCPTProcedure AS proce
  ON rad.RadiologyProcedureSID=proce.RadiologyCPTProcedureSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON rad.PatientSID=pat.PatientSID
INNER JOIN CDWWork.Dim.Location AS loc
  ON rad.RequestingLocationSID=loc.LocationSID
INNER JOIN CDWWork.Dim.StopCode AS sc
  ON loc.PrimaryStopCodeSID=sc.StopCodeSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND rad.ExamDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, rad.ExamDateTime)<=(365)
AND LocationType='CLINIC' AND NOT StopCode='130'
GROUP BY cohort.PatientICN
")

radiology <- data.table(radiology)
pcp_claims <- merge(x=pcp_claims, y=radiology, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$N_imaging[which(is.na(pcp_claims$N_imaging))] <- 0
pcp_claims$N_imaging[which(pcp_claims$N_imaging>=10)] <- 10


Ref <- felm(Ref~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref_MH <- felm(Ref_MH~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
               +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref_PH <- felm(Ref_PH~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
               +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_panel <- felm(N_panel~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_imaging <- felm(N_imaging~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                  +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

mh <- stargazer(Ref, Ref_MH, Ref_PH, N_panel, N_imaging,
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Any Referral", "MH Ref", "Cardiology Ref", "N Panels", "N Imaging"))

Ref <- felm(Ref~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref_MH <- felm(Ref_MH~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
               +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref_PH <- felm(Ref_PH~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
               +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_panel <- felm(N_panel~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_imaging <- felm(N_imaging~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                  +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

ph <- stargazer(Ref, Ref_MH, Ref_PH, 
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("Any Referral", "MH Ref", "Cardiology Ref", "Cardiology Ref", "N Panels", "N Imaging"))

Ref <- felm(Ref~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
            +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref_MH <- felm(Ref_MH~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
               +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
Ref_PH <- felm(Ref_PH~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
               +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_panel <- felm(N_panel~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
N_imaging <- felm(N_imaging~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                  +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

acsc <- stargazer(Ref, Ref_MH, Ref_PH, 
                  apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=2, dep.var.labels=c("Any Referral", "MH Ref", "Cardiology Ref", "Cardiology Ref", "N Panels", "N Imaging"),
                  add.lines=list(c("Mean Dep. Var. (x100)", 100*round(mean(Ref$fitted.values),3), 100*round(mean(Ref_MH$fitted.values),3), 100*round(mean(Ref_PH$fitted.values),3))))

rbind(mh[1:17], ph[15:17], acsc[15:24])
rm(Ref, Ref_MH, Ref_PH, referrals, N_panel, N_imaging, laboratory, radiology)




############################################################################
################                  TABLE 7               ####################
############################################################################

##### MH GUIDELINES (PANEL A)
surveys <- sqlQuery(db.handle, "set nocount on

SELECT DISTINCT
  cohort.PatientICN,
  CASE 
    WHEN SurveyName IN ('PHQ9', 'PHQ-2') THEN 'MDD'
    WHEN SurveyName IN ('PC PTSD', 'PC-PTSD', 'PCL-5', 'PCLC', 'PCLM') THEN 'PTSD'
    WHEN SurveyName IN ('AUDC') THEN 'SUD'
  END AS screen
FROM CDWWork.MH.SurveyResult AS survey
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON survey.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND survey.SurveyGivenDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, survey.SurveyGivenDateTime)<=365
WHERE SurveyGivenDateTime>=CAST('2008-10-01 00:00:00' AS datetime2(0)) AND SurveyGivenDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
AND SurveyName IN ('PHQ9', 'PHQ-2', 'PC PTSD', 'PC-PTSD', 'PCL-5', 'PCLC', 'PCLM', 'AUDC')
")

pcp_claims$MDD_screen <- (pcp_claims$PatientICN %in% surveys$PatientICN[which(surveys$screen=="MDD")]); pcp_claims$MDD_screen[which(pcp_claims$VisitDateTime<as.Date('2008-10-01'))] <- NA
pcp_claims$PTSD_screen <- (pcp_claims$PatientICN %in% surveys$PatientICN[which(surveys$screen=="PTSD")]); pcp_claims$PTSD_screen[which(pcp_claims$VisitDateTime<as.Date('2008-10-01'))] <- NA
pcp_claims$SUD_screen <- (pcp_claims$PatientICN %in% surveys$PatientICN[which(surveys$screen=="SUD")]); pcp_claims$SUD_screen[which(pcp_claims$VisitDateTime<as.Date('2008-10-01'))] <- NA


MDD_screen <- felm(MDD_screen~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
PTSD_screen <- felm(PTSD_screen~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                    +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
SUD_screen <- felm(SUD_screen~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

mh <- stargazer(MDD_screen, PTSD_screen, SUD_screen, 
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("MDD", "PTSD", "SUD"))

MDD_screen <- felm(MDD_screen~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
PTSD_screen <- felm(PTSD_screen~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                    +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
SUD_screen <- felm(SUD_screen~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

ph <- stargazer(MDD_screen, PTSD_screen, SUD_screen, 
                apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("MDD", "PTSD", "SUD"))

MDD_screen <- felm(MDD_screen~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
PTSD_screen <- felm(PTSD_screen~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                    +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
SUD_screen <- felm(SUD_screen~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount
                   +PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

acsc <- stargazer(MDD_screen, PTSD_screen, SUD_screen, 
                  apply.coef=function(x){x*100}, apply.se=function(x){x*100}, keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=2, dep.var.labels=c("MDD", "PTSD", "SUD"),
                  add.lines=list(c("Mean Dep. Var. (x100)", 100*round(mean(MDD_screen$fitted.values),3), 100*round(mean(PTSD_screen$fitted.values),3), 100*round(mean(SUD_screen$fitted.values),3))))

rbind(mh[1:17], ph[15:17], acsc[15:25])
rm(MDD_screen, PTSD_screen, SUD_screen, surveys)



##### PHYSICAL GUIDELINES (PANEL B)

# CRC: Comes from 5 different sources:
#1. Outpat CPT codes
#2. Chem lab
#3. Consult (referral) request item
#4. CPRS orderable item
#5. Radiology

colon <- sqlQuery(db.handle," set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #colon_cpt; DROP TABLE IF EXISTS #colon_lab; DROP TABLE IF EXISTS #colon_service; DROP TABLE IF EXISTS #colon_orderableitem; DROP TABLE IF EXISTS #colon_hf; DROP TABLE IF EXISTS #colon_radiology

SELECT DISTINCT CPTSID, 
CASE WHEN CPTCode IN ('G0328', 'G0464', '82270', '82274', 'G0107') OR CPTDescription LIKE '%OCCULT%' THEN 'FOBT' ELSE 'OTHER' END AS type, 
CPTCode INTO #colon_cpt 
FROM CDWWork.Dim.CPT
WHERE (CPTCode IN ('G0328', 'G0464', '82270', '82274', 'G0107') OR CPTDescription LIKE '%OCCULT%' -- fecal occult blood test(FOBT)
OR (CPTDescription LIKE '%SIGMOID%' AND CPTDescription LIKE '%FLEXIBLE%') --flexible sigmoidoscopy
OR CPTCode IN ('G0106', 'G0120', 'G0122') -- double-contrast barium enema
OR (CPTDescription LIKE '%COLONOS%')
OR ((TRY_CAST(CPTCode AS INT)>=45378 AND TRY_CAST(CPTCode AS INT)<=45398) OR CPTCode IN ('G0105', 'G0120', 'G0121', '0066T', '74263')) --colonoscopy, including virtual (colonography)
) AND CPTDescription NOT LIKE '%DIAG%'

SELECT DISTINCT a.LabChemTestSID, 'FOBT' AS type INTO #colon_lab 
FROM CDWWork.Dim.LabChemTest AS a
INNER JOIN CDWWork.Dim.NationalVALabCode AS b
  ON a.NationalVALabCodeSID=b.NationalVALabCodeSID
WHERE (LabProcedure LIKE '%OCCULT%'
OR LabProcedure LIKE '%FOBT%'
OR LabChemTestName LIKE '%FOBT%'
OR LabChemTestName LIKE '%OCCULT%'
) AND LabChemTestName NOT LIKE '%DIAG%' AND LabProcedure NOT LIKE '%DIAG%'

SELECT DISTINCT RequestServiceSID,
CASE WHEN ServiceName LIKE '%OCCULT%' OR ServiceName LIKE '%FOBT%' THEN 'FOBT' ELSE 'OTHER' END AS type INTO #colon_service
FROM CDWWork.Dim.RequestService
WHERE (ServiceName LIKE '%OCCULT%'
OR ServiceName LIKE '%FOBT%'
OR (ServiceName LIKE '%SIGMOID%' AND ServiceName LIKE '%FLEXIBLE%')
OR (ServiceName LIKE '%BARIUM%' AND ServiceName LIKE '%ENEMA%')
OR ServiceName LIKE '%COLONOS%'
OR ServiceName LIKE '%COLONOG%')
AND ServiceName NOT LIKE '%DIAG%'

SELECT DISTINCT OrderableItemSID,
CASE WHEN OrderableItemName LIKE '%OCCULT%' OR OrderableItemName LIKE '%FOBT%' THEN 'FOBT' ELSE 'OTHER' END AS type
INTO #colon_orderableitem 
FROM CDWWork.Dim.OrderableItem
WHERE (OrderableItemName LIKE '%OCCULT%'
OR OrderableItemName LIKE '%FOBT%'
OR (OrderableItemName LIKE '%SIGMOID%' AND OrderableItemName LIKE '%FLEXIBLE%')
OR (OrderableItemName LIKE '%BARIUM%' AND OrderableItemName LIKE '%ENEMA%')
OR OrderableItemName LIKE '%COLONOS%'
OR OrderableItemName LIKE '%COLONOG%')
AND OrderableItemName NOT LIKE '%DIAG%'

SELECT DISTINCT HealthFactorTypeSID, 'FOBT' AS type INTO #colon_hf 
FROM CDWWork.Dim.HealthFactorType WHERE HealthFactorType LIKE '%FOBT%'

SELECT DISTINCT RadiologyProcedureSID, 
CASE WHEN RadiologyProcedure LIKE '%OCCULT%' OR RadiologyProcedure LIKE '%FOBT%' THEN 'FOBT' ELSE 'OTHER' END AS type
INTO #colon_radiology 
FROM CDWWork.Dim.RadiologyProcedure
WHERE (RadiologyProcedure LIKE '%OCCULT%'
OR RadiologyProcedure LIKE '%FOBT%'
OR (RadiologyProcedure LIKE '%SIGMOID%' AND RadiologyProcedure LIKE '%FLEXIBLE%')
OR (RadiologyProcedure LIKE '%BARIUM%' AND RadiologyProcedure LIKE '%ENEMA%')
OR RadiologyProcedure LIKE '%COLONOS%'
OR RadiologyProcedure LIKE '%COLONOG%')
AND RadiologyProcedure NOT LIKE '%DIAG%'
UNION(
SELECT RadiologyProcedureSID, CASE WHEN #colon_cpt.CPTCode IN ('G0328', 'G0464', '82270', '82274', 'G0107') THEN 'FOBT' ELSE 'OTHER' END AS type
FROM CDWWork.Dim.RadiologyProcedure AS a
INNER JOIN #colon_cpt
  ON a.CPTSID=#colon_cpt.CPTSID
)


SELECT DISTINCT
  cohort.PatientICN, type
FROM CDWWork.Outpat.VProcedure AS visit
INNER JOIN #colon_cpt
  ON visit.CPTSID=#colon_cpt.CPTSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.VisitDateTime)<=365
WHERE visit.VisitDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND visit.VisitDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
UNION(
SELECT 
  cohort.PatientICN, type
FROM CDWWork.Chem.PatientLabChem AS chem
INNER JOIN #colon_lab
  ON chem.LabChemTestSID=#colon_lab.LabChemTestSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON chem.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND chem.LabChemSpecimenDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, chem.LabChemSpecimenDateTime)<=365
WHERE chem.LabChemSpecimenDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND chem.LabChemSpecimenDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT 
  cohort.PatientICN, type
FROM CDWWork.Con.Consult AS con
INNER JOIN #colon_service
  ON con.ToRequestServiceSID=#colon_service.RequestServiceSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON con.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND con.RequestDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, con.RequestDateTime)<=365
WHERE con.RequestDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND con.RequestDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT 
  cohort.PatientICN, type
FROM CDWWork.CPRSOrder.OrderedItem AS orderitem
INNER JOIN #colon_orderableitem
  ON orderitem.OrderableItemSID=#colon_orderableitem.OrderableItemSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON orderitem.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND orderitem.OrderStartDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, orderitem.OrderStartDateTime)<=365
WHERE orderitem.OrderStartDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND orderitem.OrderStartDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT 
  cohort.PatientICN, type
FROM CDWWork.HF.HealthFactor AS hf
INNER JOIN #colon_hf
  ON hf.HealthFactorTypeSID=#colon_hf.HealthFactorTypeSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON hf.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND hf.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, hf.VisitDateTime)<=365
WHERE hf.VisitDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND hf.VisitDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT 
  cohort.PatientICN, type
FROM CDWWork.Rad.RadiologyExam AS radiology
INNER JOIN #colon_radiology
  ON radiology.RadiologyProcedureSID=#colon_radiology.RadiologyProcedureSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON radiology.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND radiology.ExamDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, radiology.ExamDateTime)<=365
WHERE radiology.ExamDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND radiology.ExamDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
)
")

pcp_claims$colon <- (pcp_claims$PatientICN %in% colon$PatientICN)
pcp_claims$colon[which(pcp_claims$Age<50|pcp_claims$Age>=75)] <- NA
pcp_claims$colon_fobt <- (pcp_claims$PatientICN %in% colon$PatientICN[which(colon$type=="FOBT")]); pcp_claims$colon_fobt[which(pcp_claims$Age<50|pcp_claims$Age>=75)] <- NA
pcp_claims$colon_other <- (pcp_claims$PatientICN %in% colon$PatientICN[which(colon$type=="OTHER")]); pcp_claims$colon_other[which(pcp_claims$Age<50|pcp_claims$Age>=75)] <- NA


hcv <- sqlQuery(db.handle,"set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #hcv_cpt; DROP TABLE IF EXISTS #hcv_lab; DROP TABLE IF EXISTS #hcv_service; DROP TABLE IF EXISTS #hcv_orderableitem

SELECT DISTINCT CPTSID INTO #hcv_cpt FROM CDWWork.Dim.CPT WHERE CPTCode IN ('86803', '87522', '87902', '80074', '87521', '86804', '87520', 'G0472')

SELECT DISTINCT a.LabChemTestSID INTO #hcv_lab 
FROM CDWWork.Dim.LabChemTest AS a
INNER JOIN CDWWork.Dim.NationalVALabCode AS b
  ON a.NationalVALabCodeSID=b.NationalVALabCodeSID
WHERE LabProcedure LIKE '%HCV%' OR LabProcedure LIKE '%HEP-C%' OR LabProcedure LIKE '%HEPC%' OR LabProcedure LIKE '%HEP C%' OR LabProcedure LIKE '%HEPATITIS C%'
OR LabChemTestName LIKE '%HCV%' OR LabChemTestName LIKE '%HEP-C%' OR LabChemTestName LIKE '%HEPC%' OR LabChemTestName LIKE '%HEP C%' OR LabChemTestName LIKE '%HEPATITIS C%'

SELECT DISTINCT RequestServiceSID INTO #hcv_service
FROM CDWWork.Dim.RequestService
WHERE ServiceName LIKE '%HCV%' OR ServiceName LIKE '%HEP-C%' OR ServiceName LIKE '%HEPC%' OR ServiceName LIKE '%HEP C%' OR ServiceName LIKE '%HEPATITIS C%'

SELECT DISTINCT OrderableItemSID INTO #hcv_orderableitem 
FROM CDWWork.Dim.OrderableItem
WHERE OrderableItemName LIKE '%HCV%' OR OrderableItemName LIKE '%HEP-C%' OR OrderableItemName LIKE '%HEPC%' OR OrderableItemName LIKE '%HEP C%' OR OrderableItemName LIKE '%HEPATITIS C%'

SELECT DISTINCT
  cohort.PatientICN
FROM CDWWork.Outpat.VProcedure AS visit
INNER JOIN #hcv_cpt
  ON visit.CPTSID=#hcv_cpt.CPTSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.VisitDateTime)<=(365)
WHERE visit.VisitDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND visit.VisitDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
UNION(
SELECT 
  cohort.PatientICN
FROM CDWWork.Chem.PatientLabChem AS chem
INNER JOIN #hcv_lab
  ON chem.LabChemTestSID=#hcv_lab.LabChemTestSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON chem.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND chem.LabChemSpecimenDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, chem.LabChemSpecimenDateTime)<=(365)
WHERE chem.LabChemSpecimenDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND chem.LabChemSpecimenDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT 
  cohort.PatientICN
FROM CDWWork.Con.Consult AS con
INNER JOIN #hcv_service
  ON con.ToRequestServiceSID=#hcv_service.RequestServiceSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON con.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND con.RequestDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, con.RequestDateTime)<=(365)
WHERE con.RequestDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND con.RequestDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT 
  cohort.PatientICN
FROM CDWWork.CPRSOrder.OrderedItem AS orderitem
INNER JOIN #hcv_orderableitem
  ON orderitem.OrderableItemSID=#hcv_orderableitem.OrderableItemSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON orderitem.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND orderitem.OrderStartDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, orderitem.OrderStartDateTime)<=(365)
WHERE orderitem.OrderStartDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND orderitem.OrderStartDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
)
")

pcp_claims$hcv <- (pcp_claims$PatientICN %in% hcv$PatientICN); pcp_claims$hcv[which(pcp_claims$Age>=79)] <- NA


hiv <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #hiv_cpt; DROP TABLE IF EXISTS #hiv_lab; DROP TABLE IF EXISTS #hiv_service; DROP TABLE IF EXISTS #hiv_orderableitem

SELECT DISTINCT CPTSID INTO #hiv_cptFROM CDWWork.Dim.CPT WHERE CPTCode IN ('86703', '86689', '86701', '87390', '87536', '86702', '86689', 'G0432', '87391', 'S3645', 'G0433', '87525', '87901', 'G0435')

SELECT DISTINCT a.LabChemTestSID INTO #hiv_lab 
FROM CDWWork.Dim.LabChemTest AS a
INNER JOIN CDWWork.Dim.NationalVALabCode AS b
  ON a.NationalVALabCodeSID=b.NationalVALabCodeSID
WHERE LabProcedure LIKE '%HIV%' OR LabChemTestName LIKE '%HIV%'

SELECT DISTINCT RequestServiceSID INTO #hiv_service
FROM CDWWork.Dim.RequestService
WHERE ServiceName LIKE '%HIV%'

SELECT DISTINCT OrderableItemSID INTO #hiv_orderableitem 
FROM CDWWork.Dim.OrderableItem
WHERE OrderableItemName LIKE '%HIV%'

SELECT DISTINCT
  cohort.PatientICN
FROM CDWWork.Outpat.VProcedure AS visit
INNER JOIN #hiv_cpt
  ON visit.CPTSID=#hiv_cpt.CPTSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.VisitDateTime)<=(365)
WHERE visit.VisitDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND visit.VisitDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
UNION(
SELECT 
  cohort.PatientICN
FROM CDWWork.Chem.PatientLabChem AS chem
INNER JOIN #hiv_lab
  ON chem.LabChemTestSID=#hiv_lab.LabChemTestSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON chem.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND chem.LabChemSpecimenDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, chem.LabChemSpecimenDateTime)<=(365)
WHERE chem.LabChemSpecimenDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND chem.LabChemSpecimenDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT 
  cohort.PatientICN
FROM CDWWork.Con.Consult AS con
INNER JOIN #hiv_service
  ON con.ToRequestServiceSID=#hiv_service.RequestServiceSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON con.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND con.RequestDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, con.RequestDateTime)<=(365)
WHERE con.RequestDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND con.RequestDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
)
UNION(
SELECT 
  cohort.PatientICN
FROM CDWWork.CPRSOrder.OrderedItem AS orderitem
INNER JOIN #hiv_orderableitem
  ON orderitem.OrderableItemSID=#hiv_orderableitem.OrderableItemSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON orderitem.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND orderitem.OrderStartDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, orderitem.OrderStartDateTime)<=(365)
WHERE orderitem.OrderStartDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND orderitem.OrderStartDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
)
")

pcp_claims$hiv <- (pcp_claims$PatientICN %in% hiv$PatientICN); pcp_claims$hiv[which(pcp_claims$Age>=64)] <- NA


flu <- sqlQuery(db.handle," set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #flu

SELECT DISTINCT ImmunizationNameSID INTO #flu FROM CDWWork.Dim.ImmunizationName WHERE CVXCode IN ('88', '140', '141', '135', '153', '144', '158', '111', '150', '149', '155', '151')
OR ImmunizationName IN ('FLU TRIVALENT (HISTORICAL)', 'Influenza (HISTORICAL)', 'Influenza Immunization (HISTORICAL)')

SELECT DISTINCT 
 cohort.PatientICN
FROM CDWWork.Immun.Immunization AS immun
INNER JOIN #flu
  ON immun.ImmunizationNameSID=#flu.ImmunizationNameSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON immun.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND immun.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, immun.VisitDateTime)<=(375)
WHERE immun.VisitDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND immun.VisitDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
")

pcp_claims$flu <- (pcp_claims$PatientICN %in% flu$PatientICN)


tobacco <- sqlQuery(db.handle,"set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #tobacco;

SELECT DISTINCT HealthFactorTypeSID INTO #tobacco FROM CDWWork.Dim.HealthFactorType  WHERE HealthFactorType LIKE '%TOB%' OR HealthFactorType LIKE '%SMO%' OR HealthFactorCategory LIKE '%TOB%' OR HealthFactorCategory LIKE '%SMO%'

SELECT DISTINCT 
  cohort.PatientICN
FROM CDWWork.HF.HealthFactor AS hf
INNER JOIN #tobacco
  ON hf.HealthFactorTypeSID=#tobacco.HealthFactorTypeSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON hf.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND hf.HealthFactorDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, hf.HealthFactorDateTime)<=(365)
WHERE hf.HealthFactorDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND hf.HealthFactorDateTime<CAST('2018-10-01 00:00:00' AS datetime2(0))
")

tobacco <- data.table(tobacco)
pcp_claims$tobacco <- (pcp_claims$PatientICN %in% tobacco$PatientICN)



colon <- felm(colon_fobt~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
hcv <- felm(hcv~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
hiv <- felm(hiv~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
flu <- felm(flu~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
tobacco <- felm(tobacco~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

mh <- stargazer(colon, hcv, hiv, flu, tobacco, apply.coef=function(x){x*100}, apply.se=function(x){x*100},
                keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("CRC", "HCV", "HIV", "Flu", "Tobacco"))

colon <- felm(colon_fobt~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
hcv <- felm(hcv~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
hiv <- felm(hiv~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
flu <- felm(flu~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
tobacco <- felm(tobacco~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

ph <- stargazer(colon, hcv, hiv, flu, tobacco, apply.coef=function(x){x*100}, apply.se=function(x){x*100},
                keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=2, dep.var.labels=c("CRC", "HCV", "HIV", "Flu", "Tobacco"))

colon <- felm(colon_fobt~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
hcv <- felm(hcv~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
hiv <- felm(hiv~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
flu <- felm(flu~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
tobacco <- felm(tobacco~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)

acsc <- stargazer(colon, hcv, hiv, flu, tobacco, apply.coef=function(x){x*100}, apply.se=function(x){x*100},
                  keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=2, dep.var.labels=c("CRC", "HCV", "HIV", "Flu", "Tobacco"),
                  add.lines=list(c("Mean Dep. Var. (x100)", 100*round(mean(colon$fitted.values),3), 100*round(mean(hcv$fitted.values),3),
                                   100*round(mean(hiv$fitted.values),3), 100*round(mean(flu$fitted.values),3), 100*round(mean(tobacco$fitted.values),3))))

rbind(mh[1:17], ph[15:17], acsc[15:25])
rm(colon, hcv, hiv, flu, tobacco)



######################################################################
#############                  TABLE 8               ################# 
######################################################################

ValueAdd <- fread("~/ValueAdd.csv")
ValueAdd$StaffSSN <- as.numeric(as.character(ValueAdd$StaffSSN))

age_gender <- sqlQuery(db.handle, "
SELECT DISTINCT
  cohort.StaffSSN, CAST(DateOfBirth AS DATE) AS DateOfBirth, staff.Gender
FROM CDWWORK.SStaff.SStaff AS staff
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON staff.StaffSSN=cohort.StaffSSN")
age_gender <- data.table(age_gender)
age_gender$StaffSSN <- as.character(age_gender$StaffSSN)
age_gender <- age_gender[!is.na(StaffSSN)]
age_gender$DateOfBirth <- as.Date(age_gender$DateOfBirth, "%Y-%m-%d")
age_gender$DateOfBirth[which(age_gender$DateOfBirth==as.Date("1901-01-01", "%Y-%m-%d"))] <- NA
age_gender$Gender <- as.character(age_gender$Gender)
age_gender$Female <- (age_gender$Gender=="FEMALE")
age_gender$Female[which(!age_gender$Gender %in% c("FEMALE", "MALE"))] <- NA

age_gender <- age_gender[, .(Female=(sum(Female, na.rm=T)>0),
                             DateOfBirth=max(DateOfBirth, na.rm=T)), by="StaffSSN"]
age_gender$StaffSSN <- as.numeric(as.character(age_gender$StaffSSN))

ValueAdd <- merge(x=ValueAdd, y=age_gender[,c("StaffSSN", "Female", "DateOfBirth")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)
ValueAdd$Age_AtFreqYear <- ValueAdd$FrequentYear-year(ValueAdd$DateOfBirth)

N_per_year <- pcp_claims[, .(N_year=length(VisitSID)), by=c("StaffSSN", "VisitYear")]
N_per_year$StaffSSN <- as.numeric(as.character(N_per_year$StaffSSN))
N_per_year <- merge(x=N_per_year, y=age_gender[,c("StaffSSN", "DateOfBirth")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)
N_per_year$age_int <- (as.numeric(N_per_year$VisitYear)-year(N_per_year$DateOfBirth))*N_per_year$N_year
N_per_year <- N_per_year[, .(Age_avg=sum(age_int)/sum(N_year)), by="StaffSSN"]
N_per_year$StaffSSN <- as.numeric(as.character(N_per_year$StaffSSN))
ValueAdd <- merge(x=ValueAdd, y=N_per_year, by.x="StaffSSN", by.y="StaffSSN", all.x=T)

ValueAdd$PhysicianDummy <- (ValueAdd$ProviderType %in% c("Allopathic & Osteopathic Physicians", "Physicians (M.D. and D.O.)"))

##### Get FTE 
fte <- sqlQuery(odbcDriverConnect('server=vhacdwa01.vha.med.va.gov; database=CDWWork; driver={SQL server}; trusted_connection=true'),
"set nocount on;

SELECT
  StaffSSN, year, COUNT(DISTINCT VisitDate) AS days
FROM(
SELECT DATEPART(year, visit.VisitDateTime) AS year, DATEPART(dayofyear, visit.VisitDateTime) AS VisitDate, cohort.StaffSSN
FROM CDWWork.Outpat.VProvider AS visit
INNER JOIN CDWWork.SStaff.SStaff AS staff
  ON visit.ProviderSID=staff.StaffSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
ON staff.StaffSSN=cohort.StaffSSN
WHERE visit.VisitDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND visit.VisitDateTime<CAST('2017-03-01 00:00:00' AS datetime2(0))
) AS z
GROUP BY StaffSSN, year")

# 0 to 1 FTE across all years we observe them in the data
fte <- data.table(fte)
fte$StaffSSN <- as.numeric(as.character(fte$StaffSSN))
fte$fte <- fte$days/250
fte$fte[which(fte$fte>=1)] <- 1
fte$PartTime <- (fte$days<240)
fte <- fte[, .(FTE=mean(fte), PartTime=mean(PartTime)), by="StaffSSN"]
ValueAdd <- merge(x=ValueAdd, y=fte[,c("StaffSSN", "FTE", "PartTime")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)


patientdays <- sqlQuery(odbcDriverConnect('server=vhacdwa01.vha.med.va.gov; database=CDWWork; driver={SQL server}; trusted_connection=true'),
                        "set nocount on;
SELECT
  StaffSSN, COUNT(*) AS patientdays, COUNT(DISTINCT VisitDate) AS days
FROM(
SELECT DISTINCT cohort.StaffSSN, visit.PatientSID, CAST(visit.VisitDateTime AS date) AS VisitDate
FROM CDWWork.Outpat.VProvider AS visit
INNER JOIN CDWWork.SStaff.SStaff AS staff
  ON visit.ProviderSID=staff.StaffSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
ON staff.StaffSSN=cohort.StaffSSN
WHERE visit.VisitDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0)) AND visit.VisitDateTime<CAST('2017-03-01 00:00:00' AS datetime2(0))
) AS z
GROUP BY StaffSSN")
patientdays$StaffSSN <- as.numeric(as.character(patientdays$StaffSSN))
patientdays$pat_perday <- patientdays$patientdays/patientdays$days
ValueAdd <- merge(x=ValueAdd, y=patientdays[,c("StaffSSN", "pat_perday")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)


# Share of MH visits that are PCMHI
pcmhi <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;

SELECT
  PatientICN, COUNT(DISTINCT VisitDate) AS N_MH, COUNT(DISTINCT CASE WHEN StopCode IN ('534', '539') THEN VisitDate ELSE NULL END) AS N_PCMHI
FROM(
SELECT DISTINCT
  cohort.PatientICN, CAST(visit.VisitDateTime AS DATE) AS VisitDate, StopCode
FROM CDWWork.Outpat.Visit AS visit
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN OMHO_QFR.ECON.appointments AS cohort
  ON pat.PatientICN=cohort.PatientICN AND visit.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.AppointmentRequestDate, visit.VisitDateTime)<=(365)
INNER JOIN CDWWork.Dim.StopCode AS sc
  ON visit.PrimaryStopCodeSID=sc.StopCodeSID
WHERE visit.VisitDateTime>=CAST('2005-01-01' AS datetime2(0)) AND Visit.VisitDateTime<CAST('2018-03-01' AS datetime2(0)) AND StopCodeName NOT LIKE '%ADMIN%' AND StopCodeName NOT LIKE '%TELE%' AND (StopCode>=502 AND StopCode<=599)
) AS z
GROUP BY PatientICN")

pcmhi <- data.table(pcmhi)
pcmhi <- merge(x=pcmhi, y=pcp_claims[,c("PatientICN", "StaffSSN")], by.x="PatientICN", by.y="PatientICN", all.y=T)
pcmhi$N_PCMHI[which(is.na(pcmhi$N_PCMHI))] <- 0; pcmhi$N_MH[which(is.na(pcmhi$N_MH))] <- 0;
pcmhi <- pcmhi[, .(Share_PCMHI=sum(N_PCMHI)/sum(N_MH), Share_Patients=mean(N_PCMHI>0)), by="StaffSSN"]
ValueAdd <- merge(x=ValueAdd, y=pcmhi[,c("StaffSSN", "Share_PCMHI", "Share_Patients")], by.x="StaffSSN", by.y="StaffSSN", all.x=T)

ValueAdd$Age_bins <- "0"
ValueAdd$Age_bins[which(is.na(ValueAdd$Age_avg))] <- "99"
ValueAdd$Age_bins[which(ValueAdd$Age_avg<25)] <- "99"
ValueAdd$Age_bins[which(ValueAdd$Age_avg>=25 & ValueAdd$Age_avg<35)] <- "1"
ValueAdd$Age_bins[which(ValueAdd$Age_avg>=35 & ValueAdd$Age_avg<45)] <- "2"
ValueAdd$Age_bins[which(ValueAdd$Age_avg>=45 & ValueAdd$Age_avg<55)] <- "3"
ValueAdd$Age_bins[which(ValueAdd$Age_avg>=55)] <- "4"

ValueAdd$Female_cat <- "0"
ValueAdd$Female_cat[which(ValueAdd$Female==1)] <- "1"
ValueAdd$Female_cat[which(is.na(ValueAdd$Female))] <- "2"

years <- pcp_claims[,.(N_Years=length(unique(VisitYear))), by="StaffSSN"]
years$StaffSSN <- as.numeric(as.character(years$StaffSSN))
ValueAdd <- merge(x=ValueAdd, y=years, by.x="StaffSSN", by.y="StaffSSN", all.x=T)

ValueAdd$N_peryear <- ValueAdd$N/ValueAdd$N_Years

# multivariate regression
mh_cor <- felm(FE_MH_std~Female_cat+Age_bins+PartTime+Share_PCMHI+pat_perday+N_peryear+PhysicianDummy|Sta3n|0|Sta3n, data=ValueAdd, weights=ValueAdd$N)
ph_cor <- felm(FE_PH_std~Female_cat+Age_bins+PartTime+Share_PCMHI+pat_perday+N_peryear+PhysicianDummy|Sta3n|0|Sta3n, data=ValueAdd, weights=ValueAdd$N)
acsc_cor <- felm(FE_ACSC_std~Female_cat+Age_bins+PartTime+Share_PCMHI+pat_perday+N_peryear+PhysicianDummy|Sta3n|0|Sta3n, data=ValueAdd, weights=ValueAdd$N)

stargazer(mh_cor, ph_cor, acsc_cor, covariate.labels=c("Physician", "Female", "Age: 35-44", "Age: 45-54", "Age: 55+", "DELETEROW", "Part Time", "Share PCMHI", "Patients per day", "New Patients Per Year"))
rm(mh_cor, ph_cor, acsc_cor)

table(rep(ValueAdd$Female_cat[which(!ValueAdd$Female_cat=="3")], ValueAdd$N[which(!ValueAdd$Female_cat=="3")]))
table(rep(ValueAdd$Age_bins[which(!ValueAdd$Age_bins=="99")], ValueAdd$N[which(!ValueAdd$Age_bins=="99")]))
weighted.mean(ValueAdd$PartTime, ValueAdd$N, na.rm=T)
weighted.mean(ValueAdd$Share_PCMHI, ValueAdd$N)
weighted.mean(ValueAdd$pat_perday, ValueAdd$N, na.rm=T)
weighted.mean(ValueAdd$N_peryear, ValueAdd$N, na.rm=T)
weighted.mean(ValueAdd$PhysicianDummy, ValueAdd$N)


##########################################################
#############             TABLE 9         ################
##########################################################

# Truncate relationshiplength at 3 years
pcp_claims$relationshipLength[which(is.na(pcp_claims$relationshipLength))] <- as.numeric(as.Date("2020-09-01", "%Y-%m-%d")-as.Date(pcp_claims$RelationshipStartDate[which(is.na(pcp_claims$relationshipLength))]))

pcp_claims$length_truncate <- pcp_claims$relationshipLength
pcp_claims$length_truncate[which(pcp_claims$relationshipLength>=(365*3))] <- 365*3
pcp_claims$length_truncate[which(is.na(pcp_claims$relationshipLength))] <- 365*3

###### Intensive margin

intensive_margin <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;

SELECT
  PatientICN, COUNT(DISTINCT VisitDate) AS N_visits
FROM(
  SELECT
    cohort.PatientICN, CAST(visit.VisitDateTime AS DATE) AS VisitDate
  FROM CDWWork.Outpat.VProvider AS visit
  INNER JOIN CDWWork.SStaff.SStaff AS staff
    ON visit.ProviderSID=staff.StaffSID
  INNER JOIN CDWWork.SPatient.SPatient AS pat
    ON visit.PatientSID=pat.PatientSID
  INNER JOIN OMHO_QFR.ECON.appointments AS cohort
    ON staff.StaffSSN=cohort.StaffSSN AND pat.PatientICN=cohort.PatientICN
  WHERE visit.VisitDateTime>=cohort.AppointmentRequestDate AND DATEDIFF(day, cohort.RelationshipStartDate, visit.VisitDateTime)<=(365)
  AND visit.VisitDateTime>=CAST('2005-01-01 00:00:00' AS datetime2(0))
) AS a
GROUP BY PatientICN")

pcp_claims <- merge(x=pcp_claims, y=intensive_margin, by.x="PatientICN", by.y="PatientICN", all.x=T)
pcp_claims$N_visits[which(is.na(pcp_claims$N_visits))] <- 0
pcp_claims$N_visits[which(pcp_claims$N_visits>=20)] <- 20

# column 2 focuses on those alive for at least 3 years
length_truncate <- felm(length_truncate~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
length_truncate2 <- felm(length_truncate~FE_MH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

mh <- stargazer(length_truncate, length_truncate2,
                keep=c("FE_MH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=3, dep.var.labels=c("Relationship Length"))


length_truncate <- felm(length_truncate~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
length_truncate2 <- felm(length_truncate~FE_PH_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

ph <- stargazer(length_truncate, length_truncate2,
                keep=c("FE_PH_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                digits=3, dep.var.labels=c("Relationship Length"))


length_truncate <- felm(length_truncate~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims)
length_truncate2 <- felm(length_truncate~FE_ACSC_std+RaceEthnicity+factor(Married)+Medicare+Medicaid+PeriodOfService+ServiceConnectedFlag+UnemployableFlag+ExposedToAgentOrangeFlag+RadiationExposureIndicatedFlag+AnnualVACheckAmount+PriorVisit+MH_dummy_lag1_v1+PH_dummy_lag1_v1+ACSC_lag1+WaitTime|Age_bins+EnrollmentPriority+InstitutionSID+YearXMonth+DayOfWeek+TimeOut_bins|0|InstitutionSID, data=pcp_claims[DeathRelativeMonth>36])

acsc <- stargazer(length_truncate, length_truncate2,
                  keep=c("FE_ACSC_std"), omit.stat=c("rsq", "adj.rsq", "ser"), keep.stat=c("N"),
                  digits=3, dep.var.labels=c("Relationship Length"))

rbind(mh[1:17], ph[15:17], acsc[15:24])
rm(length_truncate, length_truncate2)